## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin Phi[2] = valeur du noeud Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,90,5),
omega2 = c(0.5,0.1,0.01),
#Survival data,
nu2 = 0.5,
a = 90,
b = 50,
alpha = 7,
beta = 10)
#=======================================#
t <- seq(60,120, length.out = 10) #value of times
set.seed(18031997)
data <- JLS_data(G = 8, ng = 4, time = t, fct = m, param = param, linkfct = function(t,eta,phi) phi[,2]/phi[,3])
getDim(data)
## G ng N n F.
## 8 4 32 320 3
Y <- data$obs
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1
Phi <- fct_vector(function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta) {
c(- n/(2*sigma2), -N/(2*rho2), #1, 2
- G/(2*omega2), mu, #3, 4
- G/(2*nu2), #5
a/sigma2_a, -1/(2*sigma2_a), #6, 7
b/sigma2_b , -1/(2*sigma2_b), #8, 9
alpha/sigma2_alpha, -1/(2*sigma2_alpha), #10, 11
beta ) },#12
dim = c(1,1,F., F., rep(1,8)) )$eval
zeta <- function(beta, eta, phi, gamma, a,b, alpha)
{
lbd <- function(t,g) t^{b-1} * exp(alpha*data@linkfct(t, eta[g], matrix(phi[g,], nrow = 1)))
P1 <- 1:nrow(data@survival) %>% sapply(function(i) lbd(data@survival$obs[i], data@survival$gen[i]) )
# sapply(function(i) integrate(function(t) lbd(t,data@survival$gen[i]), 0, data@survival$obs[i])$value )
P2 <- exp(beta*data@survival$U + as(rep(gamma, each = ng), 'matrix'))
beta*sum(data@survival$U) - b*a^-b * sum(P2*P1)
}
zeta.der <- function(beta, eta, phi, gamma, a, b, alpha)
{
lbd <- function(t,g) t^{b-1} * exp(alpha*data@linkfct(t, eta[g], matrix(phi[g,], nrow = 1)))
P1 <- 1:nrow(data@survival) %>% sapply(function(i) data@survival$obs[i]/b*lbd(data@survival$obs[i], data@survival$gen[i]) )
# sapply(function(i) integrate(function(t) lbd(t,data@survival$gen[i]), 0, data@survival$obs[i])$value )
P2 <- data@survival$U*exp(beta*data@survival$U + as(rep(gamma, each = ng), 'matrix'))
sum(data@survival$U) - b*a^-b * sum(P2*P1)
}
S <- fct_vector(function(eta, phi, ...) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
function(eta, ...) mean(eta^2), #2
function(phi, ...) apply(phi^2, 2, mean), #3
function(phi, ...) apply(phi, 2, mean), #4
function(gamma, ...) mean(gamma^2), #5
function(a, ...) a, #6
function(a, ...) a^2, #7
function(b, ...) b, #8
function(b, ...) b^2, #9
function(alpha, ...) alpha, #10
function(alpha, ...) alpha^2, #11
function(eta, phi, gamma, a, b, alpha) #12
uniroot(function(beta)zeta.der(beta, eta, phi, gamma, a,b,alpha),
lower = -1000, upper = 1000)$root,#, tol = 1e-3)$root,
dim = c(1,1,F., F., rep(1,8)) )
h <- fct_vector(function(a,b, ...) N*log(b*a^-b),
function(b, ...) (b-1)*sum(log(t)),
function(gamma,...) ng*sum(gamma),
function(eta, phi, ...) sum(data@fct(1, eta, apply(phi, 2, rep, ng) )))
S[12](data@eta, data@phi, data@gamma, param$a, param$b, param$alpha)
## [1] 4.708159
## attr(,"12")
## [1] 1
loglik.phi <- function(x, eta, gamma, a,b,alpha,Phi){
id <- c(1,3,4)
xp <- setoffset(x, Phi%a%4)
sum(Phi%a%1*S$eval(eta = eta, phi = xp, i = 1) + Phi%a%3 * S$eval(phi = x, i = 3) ) +
zeta(Phi%a%12, eta, xp, gamma, a, b, alpha)
}
loglik.eta <- function(x, phi, gamma, a,b,alpha, Phi){
id <- c(1,2)
sum( Phi%a%id * S$eval(eta = x, phi = phi, i = id) ) +
zeta(Phi%a%12, x, phi, gamma, a, b, alpha)
}
loglik.gamma <- function(x, eta, phi, a, b, alpha, Phi){
Phi%a%5 * sum(x^2) + ng*sum(x) + h$eval(gamma = x, i = 3) +
zeta(Phi%a%12, eta, phi, x, a, b, alpha) #Correction fait à la va vite
}
#c'est faux à refaire !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
loglik.a <- function(x, b, Phi) sum( Phi%a%6:7 * S$eval(a = x, b = b, i = 6:7) ) + h$eval(a = x, b = b, i = 1)
loglik.b <- function(x, a, Phi) sum( Phi%a%8:9 * S$eval(b = x, i = 8:9) ) + sum(h$eval(a = a, b = x, i = c(1,2)))
loglik.alpha <- function(x, eta, phi, gamma, a,b,Phi){
id <- 10:11
sum( Phi%a%id * S$eval(alpha = x, i = id) ) +
zeta(Phi%a%12, eta, phi, gamma, a, b, x) +
x*h$eval(eta = eta, phi = phi, i = 4)
}
# --- Initialisation des paramètres --- #
para <- param %>% sapply(function(x) x* runif(1, 1.1,1.4))
# para$rho2 = 0.2 ; para$omega2 <- rep(.1,3)
# --- Initialisation des chaines MC : Z_0 --- #
Z <- getLatente(data, format = 'list') #list()
Z$eta = list( matrix(rep(0, G*ng), ncol = 1) )
Z$phi = list( matrix(rep(para$mu, G), nrow = F.) %>% t )
Z$gamma = list(matrix(rep(0, G), ncol = 1))
Z$a = list(matrix(para$a))
Z$b = list(matrix(para$b))
Z$alpha = list(matrix(para$alpha))
sim <- function(niter, h, Phih, eta, phi, gamma, a, b, alpha, verbatim = F)
{
eta <- list(MH_High_Dim_future(niter, eta[[1]], sd = sd.eta, loglik.eta,
phi[[1]], gamma[[1]], a[[1]], b[[1]], alpha[[1]],
Phih, cores = ncores, verbatim = verbatim ))
phi <- list(MH_Gibbs_Sampler_future(niter, setoffset(phi[[1]],-Phih%a%4), sd.phi, loglik.phi,
eta[[1]], gamma[[1]], a[[1]], b[[1]], alpha[[1]],
Phih, cores = 1, verbatim = verbatim )) %>%
lapply(setoffset, Phih%a%4)
gamma <- list(MH_High_Dim_future(niter, gamma[[1]], sd = .03, loglik.gamma,
eta[[1]], phi[[1]], a[[1]], b[[1]], alpha[[1]],
Phih, cores = 1, verbatim = verbatim ))
a <- list(MH_High_Dim_future(niter, a[[1]], sd = .5, loglik.a, b[[1]],
Phih, cores = 1, verbatim = verbatim ))
b <- list(MH_High_Dim_future(niter, b[[1]], sd = .5, loglik.b, a[[1]],
Phih, cores = 1, verbatim = verbatim ))
alpha <- list( MH_High_Dim_future(niter, alpha[[1]], sd = 1, loglik.alpha,
eta[[1]], phi[[1]], gamma[[1]], a[[1]], b[[1]],
Phih, cores = 1, verbatim = verbatim ))
list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}
#res <- simulation_test(sim, Phi, param, 1000, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))
#names(res) %>% sapply(function(var) assign(var, res[[var]][[1]],envir = .GlobalEnv))
# oracle(maxi, S, data, a = param$a, b = param$b, alpha = param$alpha)
maxi <- function(S)
{
list(sigma2 = param$sigma2,#S%a%1,
rho2 = param$rho2,#S%a%2,
mu = param$mu,#S%a%4,
omega2 = param$omega2,#S%a%3 - (S%a%4)^2,
nu2 = S%a%5,
a = S%a%6, b = S%a%8,
alpha = S%a%10, beta = param$beta)#S%a%12 )
}
niter <- 20
MH.iter <- 10
sd.eta <- .3
sd.phi <- c(.05, .3, .05)
gg <- simulation_test(sim, Phi, param, 250, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))
res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi, eps = 1e-2,
burnin = 150, coef.burnin = 2/3, verbatim = 2)
saveRDS(res, 'saem.rds')
## [1] "SAEM execution time = 01min 59sec"
## $plot_parameter
##
## $plot_MCMC
##
## $plot_acceptation
## TableGrob (3 x 2) "arrange": 6 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]
| sigma2 | rho2 | mu.1 | mu.2 | mu.3 | omega2.1 | omega2.2 | omega2.3 | nu2 | a | b | alpha | beta | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Real value | 0.05 | 0.1 | 5 | 90 | 5 | 0.5 | 0.1 | 0.01 | 0.2899 | 90.0000 | 50.0000 | 7.0000 | 10 |
| Estimated value | 0.05 | 0.1 | 5 | 90 | 5 | 0.5 | 0.1 | 0.01 | 0.0000 | 74.0437 | 30.0973 | 2.9938 | 10 |
| Rrmse | 0.00 | 0.0 | 0 | 0 | 0 | 0.0 | 0.0 | 0.00 | 1.0000 | 0.1773 | 0.3981 | 0.5723 | 0 |